function [alphas,betas,R2s, mats_yrs] = CSregs_ver2_real(gesol, Phi_Nxz, stationary_prob_Nxz, FLAG_RUN_PARALLEL)
% This code implements the following version of Campbell-Shiller:
% hpxr(t->t+12,n) regressed on slope(t,y) = y(t,n) - y(t,1)

    if nargin<4; FLAG_RUN_PARALLEL = false; end
    
    mats_yrs = 2:5;
    
    alphas = zeros(size(mats_yrs));
    betas = zeros(size(mats_yrs));
    R2s = zeros(size(mats_yrs));
    
    Phi_Nxz_12mo = Phi_Nxz^12;
    
    if FLAG_RUN_PARALLEL
        
        stationary_prob_Nxz = stationary_prob_Nxz(:);
        
        rhs_Nxz_list = cell(size(mats_yrs));
        lhs1_Nxz_list = cell(size(mats_yrs));
        lhs2_Nxz_list = cell(size(mats_yrs));
        
        for i = 1:numel(mats_yrs)
            
            mat_yr = mats_yrs(i);
            mat_mth = mat_yr*12;

            idx_1y = find(gesol.real_bonds.maturities==12);
            idx_ny = find(gesol.real_bonds.maturities==12*mat_yr);
            idx_nless1y = find(gesol.real_bonds.maturities==12*(mat_yr-1));

            if isempty(idx_1y) || isempty(idx_ny) || isempty(idx_nless1y)
                error('Some yields are not available');
            end

            % compute coefficients for RHS: slope = y(n) - y(1)
            rhs_Nxz = -log(squeeze(gesol.real_bonds.prices(:,:,:,idx_ny)))/mat_yr ...
                + log(squeeze(gesol.real_bonds.prices(:,:,:,idx_1y)));

            % LHS = p(n-1,t+1) - p(n,t) + p(1,t)
            % LHS1 = - p(n,t) + p(1,t), LHS2 at h=12 is p(n-1,t+1)

            lhs1_Nxz = -log(squeeze(gesol.real_bonds.prices(:,:,:,idx_ny))) ...
                + log(squeeze(gesol.real_bonds.prices(:,:,:,idx_1y)));
            lhs2_Nxz = log(squeeze(gesol.real_bonds.prices(:,:,:,idx_nless1y)));
            
            rhs_Nxz_list{i} = rhs_Nxz(:);
            lhs1_Nxz_list{i} = lhs1_Nxz(:);
            lhs2_Nxz_list{i} = lhs2_Nxz(:);
            
        end
        
        parfor i = 1:numel(mats_yrs)

            % reg coeffs
            mean_rhs = sum(stationary_prob_Nxz.*rhs_Nxz_list{i});
            mean_lhs1 = sum(stationary_prob_Nxz.*lhs1_Nxz_list{i});
            mean_lhs2 = sum(stationary_prob_Nxz.*lhs2_Nxz_list{i});
            var_rhs = sum(stationary_prob_Nxz.*(rhs_Nxz_list{i}.^2)) - mean_rhs^2;
            var_lhs1 = sum(stationary_prob_Nxz.*(lhs1_Nxz_list{i}.^2)) - mean_lhs1^2;
            var_lhs2 = sum(stationary_prob_Nxz.*(lhs2_Nxz_list{i}.^2)) - mean_lhs2^2;
            cov_lhs1_lhs2 = sum(stationary_prob_Nxz.*lhs1_Nxz_list{i}.*(Phi_Nxz_12mo*lhs2_Nxz_list{i})) - mean_lhs1*mean_lhs2;
            var_lhs = var_lhs1 + var_lhs2 + 2*cov_lhs1_lhs2;
            cov_rhs_lhs1 = sum(stationary_prob_Nxz.*rhs_Nxz_list{i}.*lhs1_Nxz_list{i}) - mean_rhs*mean_lhs1;
            cov_rhs_lhs2 = sum(stationary_prob_Nxz.*rhs_Nxz_list{i}.*(Phi_Nxz_12mo*lhs2_Nxz_list{i})) - mean_rhs*mean_lhs2;

            betas(i) = (cov_rhs_lhs1 + cov_rhs_lhs2)/var_rhs;
            alphas(i) = mean_lhs1 + mean_lhs2 - betas(i)*mean_rhs;
            R2s(i) = (betas(i)^2)*var_rhs/var_lhs;

        end
        
    else
    
        for i = 1:numel(mats_yrs)

            mat_yr = mats_yrs(i);
            mat_mth = mat_yr*12;

            idx_1y = find(gesol.real_bonds.maturities==12);
            idx_ny = find(gesol.real_bonds.maturities==12*mat_yr);
            idx_nless1y = find(gesol.real_bonds.maturities==12*(mat_yr-1));

            if isempty(idx_1y) || isempty(idx_ny) || isempty(idx_nless1y)
                error('Some yields are not available');
            end

            % compute coefficients for RHS: slope = y(n) - y(1)
            rhs_Nxz = -log(squeeze(gesol.real_bonds.prices(:,:,:,idx_ny)))/mat_yr ...
                + log(squeeze(gesol.real_bonds.prices(:,:,:,idx_1y)));

            % LHS = p(n-1,t+1) - p(n,t) + p(1,t)
            % LHS1 = - p(n,t) + p(1,t), LHS2 at h=12 is p(n-1,t+1)

            lhs1_Nxz = -log(squeeze(gesol.real_bonds.prices(:,:,:,idx_ny))) ...
                + log(squeeze(gesol.real_bonds.prices(:,:,:,idx_1y)));
            lhs2_Nxz = log(squeeze(gesol.real_bonds.prices(:,:,:,idx_nless1y)));

            % reg coeffs
            mean_rhs = sum(stationary_prob_Nxz(:).*rhs_Nxz(:));
            mean_lhs1 = sum(stationary_prob_Nxz(:).*lhs1_Nxz(:));
            mean_lhs2 = sum(stationary_prob_Nxz(:).*lhs2_Nxz(:));
            var_rhs = sum(stationary_prob_Nxz(:).*(rhs_Nxz(:).^2)) - mean_rhs^2;
            var_lhs1 = sum(stationary_prob_Nxz(:).*(lhs1_Nxz(:).^2)) - mean_lhs1^2;
            var_lhs2 = sum(stationary_prob_Nxz(:).*(lhs2_Nxz(:).^2)) - mean_lhs2^2;
            cov_lhs1_lhs2 = sum(stationary_prob_Nxz(:).*lhs1_Nxz(:).*(Phi_Nxz_12mo*lhs2_Nxz(:))) - mean_lhs1*mean_lhs2;
            var_lhs = var_lhs1 + var_lhs2 + 2*cov_lhs1_lhs2;
            cov_rhs_lhs1 = sum(stationary_prob_Nxz(:).*rhs_Nxz(:).*lhs1_Nxz(:)) - mean_rhs*mean_lhs1;
            cov_rhs_lhs2 = sum(stationary_prob_Nxz(:).*rhs_Nxz(:).*(Phi_Nxz_12mo*lhs2_Nxz(:))) - mean_rhs*mean_lhs2;

            betas(i) = (cov_rhs_lhs1 + cov_rhs_lhs2)/var_rhs;
            alphas(i) = mean_lhs1 + mean_lhs2 - betas(i)*mean_rhs;
            R2s(i) = (betas(i)^2)*var_rhs/var_lhs;

        end
    
    end

end